Example 1: Circular scalings in 2d

<- back to main page

4.1. Example 1: Circular scalings in 2d#

As our first example we consider a two-dimensional example, namely \(\Omega:=\mathbb R^2\setminus B_{R}(\mathbf x_0)\), \(\Gamma:=\partial B_1(0)\), for some center \(\mathbf x_0\in B_1(0)\) and radius \(R\) such that \(R<1-\|\mathbf x_0\|\), and polar coordinates (cf. Remark 3.3)

(4.1)#\[\begin{align}\xi(\mathbf x)&:=\|\mathbf x\|-1,& \hat x(\mathbf x)&:=\frac{\mathbf x}{\|\mathbf x\|}.\end{align}\]

We create the according mesh and set some parameters.

from ngsolve import *
from netgen.geom2d import *
from nonlin_arnoldis import *
from ngsolve.webgui import Draw
from numpy import array,sqrt,loadtxt
from matplotlib.pyplot import plot,show,xlim,ylim,legend


N = 25         #infinite elements
maxh = 0.1     #mesh-size
sigma = 0.3+0.5j   #complex scaling paramter
order = 5      #fem order
shift = 5-0.5j      #shift for Arnoldi algorithm
center = (0.2,0)    #center of inner circle
R = 0.5            #radius of inner circle

#create geometry
geo = SplineGeometry()
geo.AddCircle((0,0),1,leftdomain=1,rightdomain=0,bc='Gamma')
geo.AddCircle(center,R,leftdomain=0,rightdomain=1,bc='inner')


#create mesh
mesh = Mesh(geo.GenerateMesh(maxh=maxh))
mesh.Curve(2*order)
Draw(mesh)
BaseWebGuiScene

The weak formulation after a re-scaling of the solution and testfunction in the exterior is given by (7.11).

For this first simple example we choose a frequency independent scaling $\(\sigma(\omega):=\sigma_0\in\mathbb C.\)$

We start by creating the large finite element space for implementing our tensor-product method.

Gamma = mesh.Boundaries('Gamma')

fes_int = H1(mesh,order=order,complex=True)
fes_surf = H1(mesh,order=order,complex=True,definedon=Gamma)

fes = ProductSpace(fes_int,*( N*[fes_surf]) )

We import the necessary infinte element matrices and prepare the radial matrices:

from infinite_elements import *

ie_mass,ie_laplace,__,ie_mass_x,ie_mass_xx,ie_laplace_x,ie_laplace_xx,_ = ie_matrices(N)

S_ie_1 = 1/sigma*ie_laplace+2*ie_laplace_x+sigma*ie_laplace_xx-sigma/4*ie_mass
S_ie_1[0,0]-=1/2

S_ie_2 = sigma*ie_mass

M_ie = sigma * (ie_mass+2*sigma*ie_mass_x+sigma**2*ie_mass_xx)

Now we can assemble our bilinear forms.

ds_g = ds(definedon=Gamma)
p,q = fes.TnT()
p_int,q_int = p[0],q[0]
S = BilinearForm(
    grad(p_int)*grad(q_int)*dx
    +sum(S_ie_1[i,j]*p[j]*q[i]*ds_g
       for i in range(N+1) for j in range(N+1) if abs(S_ie_1[i,j])>0)
    +sum(S_ie_2[i,j]*p[j].Trace().Deriv()*q[i].Trace().Deriv()*ds_g
       for i in range(N+1) for j in range(N+1) if abs(S_ie_2[i,j])>0)   
    ,symmetric=True).Assemble()

M = BilinearForm(
    -p_int*q_int*dx
    -sum(M_ie[i,j]*p[j]*q[i]*ds_g
       for i in range(N+1) for j in range(N+1) if abs(M_ie[i,j])>0)
    ,symmetric=True).Assemble()

Finally, we solve the resulting eigenvalue problem.

gf = GridFunction(fes,multidim=80)

#lam = sqrt(array(ArnoldiSolver(S.mat,M.mat,freedofs=fes.FreeDofs(),vecs=gf.vecs,shift=shift**2)))
lam = sqrt(array(PolyArnoldiSolver([S.mat,M.mat],shift**2,200,nevals=80,vecs=gf.vecs,inversetype='sparsecholesky',freedofs=fes.FreeDofs())))


plot(lam.real,lam.imag,'x',label='resonances')
plot([0,5*(1/sigma).real],[0,5*(1/sigma).imag],label='essential spectrum')

#load reference resonances from file
loaded=loadtxt('dhankel_1_zeros.out')
ref=(loaded[:,0]+1j*loaded[:,1])/R

plot(ref.real,ref.imag,'.k',label='reference')

xlim((0,10))
ylim((-5,0))
legend()
Draw(gf.components[0])
1/200
2/200
3/200
4/200
5/200
6/200
7/200
8/200
9/200
10/200
11/200
12/200
13/200
14/200
15/200
16/200
17/200
18/200
19/200
20/200
21/200
22/200
23/200
24/200
25/200
26/200
27/200
28/200
29/200
30/200
31/200
32/200
33/200
34/200
35/200
36/200
37/200
38/200
39/200
40/200
41/200
42/200
43/200
44/200
45/200
46/200
47/200
48/200
49/200
50/200
51/200
52/200
53/200
54/200
55/200
56/200
57/200
58/200
59/200
60/200
61/200
62/200
63/200
64/200
65/200
66/200
67/200
68/200
69/200
70/200
71/200
72/200
73/200
74/200
75/200
76/200
77/200
78/200
79/200
80/200
81/200
82/200
83/200
84/200
85/200
86/200
87/200
88/200
89/200
90/200
91/200
92/200
93/200
94/200
95/200
96/200
97/200
98/200
99/200
100/200
101/200
102/200
103/200
104/200
105/200
106/200
107/200
108/200
109/200
110/200
111/200
112/200
113/200
114/200
115/200
116/200
117/200
118/200
119/200
120/200
121/200
122/200
123/200
124/200
125/200
126/200
127/200
128/200
129/200
130/200
131/200
132/200
133/200
134/200
135/200
136/200
137/200
138/200
139/200
140/200
141/200
142/200
143/200
144/200
145/200
146/200
147/200
148/200
149/200
150/200
151/200
152/200
153/200
154/200
155/200
156/200
157/200
158/200
159/200
160/200
161/200
162/200
163/200
164/200
165/200
166/200
167/200
168/200
169/200
170/200
171/200
172/200
173/200
174/200
175/200
176/200
177/200
178/200
179/200
180/200
181/200
182/200
183/200
184/200
185/200
186/200
187/200
188/200
189/200
190/200
191/200
192/200
193/200
194/200
195/200
196/200
197/200
198/200
199/200
200/200
BaseWebGuiScene
_images/6e1d3b06002894209db051e0ae176c873c1d061757b1b7340ec29ab83b51553d.png

<- back to main page